#include<iostream>
#include<stdio.h>
#include<time.h>
#include<omp.h>
#include<arm_neon.h>
using namespace std;
int n=16;
float ** P;

#define NUM_THREADS 8

void Serial(float** A)
{
    for(int k=0;k<n;k++)
    {
        for(int j=k+1;j<n;j++)
        {
            A[k][j]/=A[k][k];
        }
        A[k][k]=1.0;
        for(int i=k+1;i<n;i++)
        {
            for(int j=k+1;j<n;j++)
            {
                A[i][j]-=A[i][k]*A[k][j];
            }
            A[i][k]=0.0;
        }
    }
}

void reset(float** A,float** B)
{
    for(int i=0;i<n;i++)
    {
        for(int j=0;j<n;j++)
        {
            B[i][j]=A[i][j];
        }
    }
}

void OMP(float** A)
{
    int i,j,k;
    float tmp;
    #pragma omp parallel num_threads(NUM_THREADS), private(i, j, k, tmp)
    for(k=0;k<n;k++)
    {
        #pragma omp single
        {
            tmp = A[k][k];
            for(j=k+1;j<n;j++)
            {
                A[k][j]/= tmp;
            }
            A[k][k]=1.0;
        }
        #pragma omp for
        for(i=k+1;i<n;i++)
        {
            tmp = A[i][k];
            for(j=k+1;j<n;j++)
            {
                A[i][j]-=tmp*A[k][j];
            }
            A[i][k]=0.0;
        }
    }
}
void OMP_static(float** A)
{
    int i,j,k;
    float tmp;
    #pragma omp parallel num_threads(NUM_THREADS), private(i, j, k, tmp)
    for(k=0;k<n;k++)
    {
        #pragma omp single
        {
            tmp = A[k][k];
            for(j=k+1;j<n;j++)
            {
                A[k][j]/= tmp;
            }
            A[k][k]=1.0;
        }
        #pragma omp for schedule(static,8)
        for(i=k+1;i<n;i++)
        {
            tmp = A[i][k];
            for(j=k+1;j<n;j++)
            {
                A[i][j]-=tmp*A[k][j];
            }
            A[i][k]=0.0;
        }
    }
}
void OMP_dynamic(float** A)
{
    int i,j,k;
    float tmp;
    #pragma omp parallel num_threads(NUM_THREADS), private(i, j, k, tmp)
    for(k=0;k<n;k++)
    {
        #pragma omp single
        {
            tmp = A[k][k];
            for(j=k+1;j<n;j++)
            {
                A[k][j]/= tmp;
            }
            A[k][k]=1.0;
        }
        #pragma omp for schedule(dynamic,8)
        for(i=k+1;i<n;i++)
        {
            tmp = A[i][k];
            for(j=k+1;j<n;j++)
            {
                A[i][j]-=tmp*A[k][j];
            }
            A[i][k]=0.0;
        }
    }
}
void OMP_guided(float** A)
{
    int i,j,k;
    float tmp;
    #pragma omp parallel num_threads(NUM_THREADS), private(i, j, k, tmp)
    for(k=0;k<n;k++)
    {
        #pragma omp single
        {
            tmp = A[k][k];
            for(j=k+1;j<n;j++)
            {
                A[k][j]/= tmp;
            }
            A[k][k]=1.0;
        }
        #pragma omp for schedule(guided,8)
        for(i=k+1;i<n;i++)
        {
            tmp = A[i][k];
            for(j=k+1;j<n;j++)
            {
                A[i][j]-=tmp*A[k][j];
            }
            A[i][k]=0.0;
        }
    }
}
void OMP_neon(float** A)
{
    int i,j,k;
    float32x4_t vt;
    float tmp;
    #pragma omp parallel num_threads(NUM_THREADS), private(i, j, k, vt,tmp)
    for(k=0;k<n;k++)
    {
        #pragma omp single
        {
            vt=vmovq_n_f32(A[k][k]);
            tmp = A[k][k];
            for(j=k+1;j+4<=n;j+=4)
            {
                float32x4_t va=vld1q_f32(&A[k][j]);
                va=vdivq_f32(va,vt);
                vst1q_f32(&A[k][j],va);
            }
            for(j;j<n;j++)
            {
                A[k][j]/=tmp;
            }
            A[k][k]=1.0;
        }
        #pragma omp for
        for(int i=k+1;i<n;i++)
        {
            float32x4_t vaik=vmovq_n_f32(A[i][k]);
            tmp = A[i][k];
            j=k+1;
            while(j%4!=0)
            {
                A[i][j]-=A[k][j]*tmp;
                j++;
            }
            for(;j+4<=n;j+=4)
            {
                float32x4_t vakj=vld1q_f32(&A[k][j]);
                float32x4_t vaij=vld1q_f32(&A[i][j]);
                float32x4_t vx=vmulq_f32(vakj,vaik);
                vaij=vsubq_f32(vaij,vx);
                vst1q_f32(&A[i][j],vaij);

            }
            for(j;j<n;j++)
            {
                A[i][j]-=A[k][j]*tmp;
            }
            A[i][k]=0.0;
        }
    }
}
void Run()
{
    float **Gauss_arr=new float*[n];
    for(int i=0;i<n;i++)
        Gauss_arr[i]=new float[n];

    for(int i=0;i<n;i++)
    {
        for(int j=0;j<i;j++)
            Gauss_arr[i][j]=0;
        Gauss_arr[i][i]=1.0;
        for(int j=i;j<n;j++)
            Gauss_arr[i][j]=rand();
    }
    for(int k=0;k<n;k++)
        for(int i=k+1;i<n;i++)
            for(int j=0;j<n;j++)
                Gauss_arr[i][j]+=Gauss_arr[k][j];

    P = new float*[n];
    for(int i=0;i<n;i++)
        P[i]=new float[n];

    struct timespec sts,ets;
    time_t dsec;
    long dnsec;
    reset(Gauss_arr,P);
    timespec_get(&sts,TIME_UTC);
    Serial(P);
    timespec_get(&ets,TIME_UTC);
    dsec=ets.tv_sec-sts.tv_sec;
    dnsec=ets.tv_nsec-sts.tv_nsec;
    if (dnsec<0)
    {
        dsec--;
        dnsec+=1000000000ll;
    }
    cout<<"Serial n: "<<n<<endl;
    printf ("%llu.%09llus\n",dsec,dnsec);
    cout<<endl;

    reset(Gauss_arr,P);
    timespec_get(&sts,TIME_UTC);
    OMP(P);
    timespec_get(&ets,TIME_UTC);
    dsec=ets.tv_sec-sts.tv_sec;
    dnsec=ets.tv_nsec-sts.tv_nsec;
    if (dnsec<0)
    {
        dsec--;
        dnsec+=1000000000ll;
    }
    cout<<"omp n: "<<n<<endl;
    printf ("%llu.%09llus\n",dsec,dnsec);
    cout<<endl;

    reset(Gauss_arr,P);
    timespec_get(&sts,TIME_UTC);
    OMP_static(P);
    timespec_get(&ets,TIME_UTC);
    dsec=ets.tv_sec-sts.tv_sec;
    dnsec=ets.tv_nsec-sts.tv_nsec;
    if (dnsec<0)
    {
        dsec--;
        dnsec+=1000000000ll;
    }
    cout<<"omp_static n: "<<n<<endl;
    printf ("%llu.%09llus\n",dsec,dnsec);
    cout<<endl;

    reset(Gauss_arr,P);
    timespec_get(&sts,TIME_UTC);
    OMP_dynamic(P);
    timespec_get(&ets,TIME_UTC);
    dsec=ets.tv_sec-sts.tv_sec;
    dnsec=ets.tv_nsec-sts.tv_nsec;
    if (dnsec<0)
    {
        dsec--;
        dnsec+=1000000000ll;
    }
    cout<<"omp_dynamic n: "<<n<<endl;
    printf ("%llu.%09llus\n",dsec,dnsec);
    cout<<endl;

    reset(Gauss_arr,P);
    timespec_get(&sts,TIME_UTC);
    OMP_guided(P);
    timespec_get(&ets,TIME_UTC);
    dsec=ets.tv_sec-sts.tv_sec;
    dnsec=ets.tv_nsec-sts.tv_nsec;
    if (dnsec<0)
    {
        dsec--;
        dnsec+=1000000000ll;
    }
    cout<<"omp_guided n: "<<n<<endl;
    printf ("%llu.%09llus\n",dsec,dnsec);
    cout<<endl;

    reset(Gauss_arr,P);
    timespec_get(&sts,TIME_UTC);
    OMP_neon(P);
    timespec_get(&ets,TIME_UTC);
    dsec=ets.tv_sec-sts.tv_sec;
    dnsec=ets.tv_nsec-sts.tv_nsec;
    if (dnsec<0)
    {
        dsec--;
        dnsec+=1000000000ll;
    }
    cout<<"omp_neon n: "<<n<<endl;
    printf ("%llu.%09llus\n",dsec,dnsec);
    cout<<endl;
}

int main()
{
    for(;n<4096;n*=2)
    {
        Run();
    }
}
